OISST Sea Surface Temperature Report
Adam A. Kemberling
Updated on: 2020-12-09
Northwest Atlantic Sea Surface Temperature Outlook
Region Extent
The spatial extent for Northwest Atlantic is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.
# File paths for various extents based on params$region
poly_path <- switch(
tolower(params$region),
"gulf of maine" = "Shapefiles/gmri_sst_focal_areas/apershing_gulf_of_maine.geojson",
"cpr gulf of maine" = "Shapefiles/gmri_sst_focal_areas/cpr_gulf_of_maine.geojson",
"northwest atlantic" = "Shapefiles/gmri_sst_focal_areas/aak_northwest_atlantic.geojson"
)
poly_path <- paste0(res_path, poly_path)
# Load the bounding box for Andy's GOM to show they align
region_extent <- st_read(poly_path, quiet = TRUE)
# Pull extents for the region for crop
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]
# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-70.875, -65.375)
crop_y <- c(40.375, 45.125)}
# Full plot
ggplot() +
#geom_stars(data = andy_stars) + # see GMRI_template_focus_areas.R
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
geom_sf(data = region_extent, color = "royalblue",
fill = "transparent", linetype = 2, size = 1) +
scale_fill_gradient2("Sea Surface Temperature Anomaly",
low = "blue",
mid = "white",
high = "red") +
map_theme +
guides("fill" = guide_colorbar(title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = T) 
Accessing OISST on Box
For many of our most-used areas climatologies and anomalies have been processed using a handful of jupyter notebook workflows to take advantage of {xarray}. The following sections will highlight what those workflows are, what they detail, and how they can be updated using the Gulf of Maine as an example.
The root folder for the following OISST products is ~/Box/NSF OKN Demo Data/oisst.
Masked Regional Timeseries
A region-masked timeseries is the most basic building block for relaying information for any of these areas. For any desired area (represented by a spatial polygon) a timeseries table of observed sea surface temperature is returned that contains the mean sea surface temperature over that area at each time step. Rather than repeat this process again for the climatology table and the anomalies these three values and a few additional fields are stored in the same table.
Pre-processed tables have been placed on box as part of the NSF C-ACCEL project under the folder: ~/Box/NSF OKN Demo Data/oisst/regional_timeseries
The naming convention in the folder is: regional_timeseries/nmfs_trawl_regions/OISST_v2_anom_ + region_name + .csv
# Switch path based on local vs. docker
timeseries_path <- switch(
EXPR = tolower(params$region),
"gulf of maine" = "oisst/regional_timeseries/gmri_sst_focal_areas/OISSTv2_anom_apershing_gulf_of_maine.csv",
"cpr gulf of maine" = "oisst/regional_timeseries/gmri_sst_focal_areas/OISSTv2_anom_cpr_gulf_of_maine.csv",
"northwest atlantic" = "oisst/regional_timeseries/gmri_sst_focal_areas/OISSTv2_anom_northwest_atlantic.csv")
timeseries_path <- str_c(okn_path, timeseries_path)
# Load timeseries
region_timeseries <- read_csv(timeseries_path,
guess_max = 1e5,
col_types = cols())
head(region_timeseries) %>% mutate_if(is_numeric, round, 2) %>% kable() %>% kable_styling()| time | modified_ordinal_day | sst | sst_clim | sst_anom | clim_sd | log_lik |
|---|---|---|---|---|---|---|
| 1981-09-01 | 244 | 10.79 | 11.16 | -0.37 | 6.49 | 2.79 |
| 1981-09-02 | 245 | 10.74 | 11.10 | -0.37 | 6.45 | 2.79 |
| 1981-09-03 | 246 | 10.76 | 11.08 | -0.32 | 6.42 | 2.78 |
| 1981-09-04 | 247 | 10.76 | 11.05 | -0.29 | 6.41 | 2.78 |
| 1981-09-05 | 248 | 10.71 | 11.01 | -0.30 | 6.41 | 2.78 |
| 1981-09-06 | 249 | 10.69 | 10.95 | -0.26 | 6.41 | 2.78 |
Climatologies are currently set up to calculate daily averages on a modified julian day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not. In these tables sst is the mean temperature observed for that date averaged across all cells within the area. sst_clim & clim_sd are the climate means and standard deviations for a 1982-2011 climatology. sst_anom is the daily observed minus the climate mean, and log_lik is the log likelihood of getting that value given a normal distribution with mean of sst_clim and sd of clim_sd.
Marine Heatwaves
The {heatwaveR} package provides a relatively quick way of working with tabular data to calculate a seasonal climate mean, as well as heatwaves and coldwaves at a desired threshold. These functions can be used to get the climatology using the standard day of year if desired, and track the number and length of heatwave events. Heatwave events follow the definition of Hobday et al. 2016
The heatwave history for the Gulf of Maine trawl region displayed above is as follows:
# Wrapper function to do heatwaves and coldwaves simultaneously at 90%
pull_heatwave_events <- function(temperature_timeseries) {
# Pull the two column dataframe for mhw estimation
test_ts <- data.frame(t = temperature_timeseries$time,
temp = temperature_timeseries$sst)
# Detect the events in a time series
ts <- ts2clm(data = test_ts, climatologyPeriod = c("1982-01-01", "2011-12-31"))
#heatwaves
mhw <- detect_event(ts)
# prep heatwave data
mhw_out <- mhw$climatology %>%
mutate(sst_anom = temp - seas) %>%
select(time = t,
sst = temp,
seas,
sst_anom,
mhw_thresh = thresh,
mhw_event = event,
mhw_event_no = event_no)
# 2. Detect cold spells
ts <- ts2clm(data = test_ts,
climatologyPeriod = c("1982-01-01", "2011-12-31"),
pctile = 10)
mcs <- detect_event(ts, coldSpells = TRUE) #coldSpells = TRUE flips boolean to < thresh
# prep cold spell data
mcs_out <- mcs$climatology %>%
select(time = t,
mcs_thresh = thresh,
mcs_event = event,
mcs_event_no = event_no)
# join heatwaves to coldwaves
hot_and_cold <- left_join(mhw_out, mcs_out, by = "time")
# 3. Data formatting for plotting,
# adds columns to plot hw and cs seperately
events_out <- hot_and_cold %>%
mutate(status = ifelse(mhw_event == TRUE, "Marine Heatwave Event", "Sea Surface Temperature"),
status = ifelse(mcs_event == TRUE, "Marine Cold Spell Event", status),
hwe = ifelse(mhw_event == TRUE, sst, NA),
cse = ifelse(mcs_event == TRUE, sst, NA),
nonevent = ifelse(mhw_event == FALSE & mcs_event == FALSE, sst, NA))
# Close the gaps between a mhw event and sst (might not need if full line for temp exists)
events_out <- events_out %>%
mutate(hwe = ifelse(is.na(hwe) & is.na(lag(hwe)) == FALSE, sst, hwe),
cse = ifelse(is.na(cse) & is.na(lag(cse)) == FALSE, sst, cse))
return(events_out)
}
# Use function
#
region_heatwaves <- pull_heatwave_events(region_timeseries)Heatwave timelines can be interactively plotted using {ggplot2} for static plots or {plotly} for interactive content.
Interactive SST Timeline
For anything we wish to host on the web there is an option to display tables and graphs that are interactive. The {plotly} package is one-such tool for producing plots that allow users to pan, zoom, and highlight discrete observations.
# Helper function
plotly_mhw_plots <- function(data){
# How to get rgb() from colors
plot_cols <- list(
"gray20" = col2rgb(col = "gray20"),
"gray40" = col2rgb(col = "gray40"),
"royalblue" = col2rgb(col = "royalblue"),
"darkred" = col2rgb(col = "darkred"),
"royalblue" = col2rgb(col = "royalblue"))
# Building the plot
fig <- plot_ly(data, x = ~time)
# Sea Surface Temperature
fig <- fig %>% add_trace(y = ~sst,
name = 'Sea Surface Temperature',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(65, 105, 225)",
width = 2))
# Heatwave Threshold
fig <- fig %>% add_trace(y = ~mhw_thresh,
name = 'MHW Threshold ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(51, 51, 51)",
width = 1,
dash = 'dot'))
# Seasonal Climatology
fig <- fig %>% add_trace(y = ~seas,
name = 'Seasonal Climatology ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(102, 102, 102)",
width = 3,
dash = 'dash'))
# Marine Cold Spell Threshold
fig <- fig %>% add_trace(y = ~mcs_thresh,
name = 'MCS Threshold ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(51, 51, 51)",
width = 1,
dash = 'dot'))
# Heatwave Event
fig <- fig %>% add_trace(y = ~hwe,
name = 'Marine Heatwave Event ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(139, 0, 0)",
width = 2))
# Cold Spell Event
fig <- fig %>% add_trace(y = ~cse,
name = 'Marine Cold Spell Event ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(65, 105, 225)",
width = 2))
# Axis Formatting
fig <- fig %>% layout(xaxis = list(title = ""),
yaxis = list (title = "Temperature (degrees C)"))
# Legend formatting
fig <- fig %>% layout(legend = list(orientation = 'h'))
return(fig)
}
# Plot the last year
last_year <- Sys.Date() - 365
region_heatwaves %>%
filter(time >= last_year) %>%
plotly_mhw_plots()Static SST Timeline
In most cases a static image will suffice, in which case {ggplot2} is a familiar go-to.
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"Cold Spell Event" = "lightblue",
"MHW Threshold" = "gray30",
"MCS Threshold" = "gray30",
"Seasonal Climatology" = "gray30")
# Set the label with degree symbol
ylab <- expression("Sea Surface Temperature"~degree~C)
# Plot the last 365 days
region_heatwaves %>%
filter(time >= last_year) %>%
ggplot(aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst), color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = cse, color = "Cold Spell Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 2, size = .5) +
geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 2, size = .5) +
geom_line(aes(y = seas, color = "Seasonal Climatology"), size = .5) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
theme(legend.title = element_blank(),
legend.position = "bottom") +
labs(x = "", y = ylab)
Warming Rates
Regional warming trends below were calculated using annual averages, and consequently the model parameters reflect annual increases in sea surface temperature.
NOTE: The trends produced below use OISST data clipped from the global extent using the “NELME” shapefile that includes Georges Bank and the Bay of Fundy.
Things to add:
* Global Trends (grey)
* Regional Trend (blue) * Coefficients on Annual scale
Annual
# Summarise by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>%
mutate(year = year(time)) %>%
group_by(year) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep")
# Plot the warming rate for the region
annual_summary %>%
ggplot(aes(year, sst_anom)) +
geom_line(group = 1) +
geom_point() +
geom_smooth(method = "lm",
aes(color = "Gulf of Maine Trend"),
formula = y ~ x, se = F, linetype = 2) +
stat_poly_eq(formula = y ~ x,
aes(color = "Gulf of Maine Trend",
label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = T) +
scale_color_manual(values = c("Gulf of Maine Trend" = "royalblue")) +
labs(x = "", y = "Sea Surface Temperature Anomaly") +
theme(legend.title = element_blank(),
legend.position = c(0.85, 0.2),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"))
Yearly Quarters
# Repeat for the seasons (equal quarters here*)
quarter_summary <- region_timeseries %>%
mutate(year = year(time),
season = factor(quarter(time, fiscal_start = 1)),
season = fct_recode(season, c("Q1" = "1"), c("Q2" = "2"), c("Q3" = "3"), c("Q4" = "4"))) %>%
group_by(year, season) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep")
quarter_summary %>%
ggplot(aes(year, sst_anom)) +
geom_line(group = 1) +
geom_point() +
geom_smooth(method = "lm",
aes(color = "Gulf of Maine Trend"),
formula = y ~ x, se = F, linetype = 2) +
stat_poly_eq(formula = y ~ x,
aes(color = "Gulf of Maine Trend",
label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = T) +
scale_color_manual(values = c("Gulf of Maine Trend" = "royalblue")) +
labs(x = "", y = "Sea Surface Temperature Anomaly") +
theme(legend.title = element_blank(),
legend.position = c(0.875, 0.05),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent")) +
facet_wrap(~season)
Heatwave Event History
As mentioned above heatwave events were determined using the methods of Hobday et al. 2016 and implemented using {heatwaveR} with a 90% threshold. The region used is again that same “NELME” Gulf of Maine region, the bottom right panel of the maps at the beginning of the document.
# Prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)
# get new axis dimensions, y = year, x = day within year
# use flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# assemble plot
ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2020.5), expand = c(0,0)) +
labs(x = "", y = "") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
Anomaly Maps
For demonstration purposes I will only be loading the 2020 global anomaly data for the below images.
# Access information to netcdf on box
okn_path <- shared.path(group = "NSF OKN", folder = "oisst/annual_anomalies/")
nc_year <- "2020"
anom_path <- str_c(okn_path, "daily_anoms_", nc_year, ".nc")
# Load 2020 as stack
anoms_2020 <- stack(anom_path)Mean Global
The following code will subset the anomalies for July and plot the average sea surface temperature anomalies for that month:
# Get the mean temperature anomalies for July
july_dates <- which(str_sub(names(anoms_2020), 7, 8) == "07")
july_avg <- mean(anoms_2020[[july_dates]])
# Convert wgs84 raster to stars array
july_st <- st_as_stars(rotate(july_avg))
# Plot global map
ggplot() +
geom_stars(data = july_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_gradient2(low = "blue",
mid = "white",
high = "red") +
map_theme +
coord_sf(expand = FALSE) +
guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Regional Map
Using the July average sst anomaly we can then clip down to the Gulf of Maine for an aerial view of the region. For this execution the bounding box of lat-lon coordinates from the regional extent will be used.
# Clip Raster - Convert to stars
shape_extent <- c(crop_x, crop_y)
region_ras <- crop(rotate(july_avg), extent(shape_extent))
region_st <- st_as_stars(region_ras)
# Get crop bounds for coord_sf
crop_x <- st_bbox(region_st)[c(1,3)]
crop_y <- st_bbox(region_st)[c(2,4)]
# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-71, -65.5)
crop_y <- c(40.5, 45)}
# Plot Gulf of Maine - wgs84
ggplot() +
geom_stars(data = region_st) +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red", na.value = "transparent") +
map_theme +
coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
guides("fill" = guide_colorbar(
title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Regionally Projected Map
These lat/lon region extent boxes are good showcases of how something that is “rectangular” by its lat/lon dimensions that has its shape change to when projected onto a non-flat surface of the earth.
To get the OISST anomaly grid to display on a curviliniear projection the originals were cropped using the corner coordinates of the extent plotted above. The raster was then warped into a regular grid using an appropriate projection as displayed below.
For the most dramatic case set the region parameter to “Northwest Atlantic”
#### Candidate coordinate reference systems ####
#### NOTE:
# transformation from rasters different than polygons
# as cells need to stretch and bend ####
# Resource: Tranformation vs. Warping
# https://r-spatial.github.io/stars/articles/stars5.html
# Robinson projection
robinson_proj <- "+proj=robin"
# Albers equal area: centered on -70 degrees
# The settings for lat_1 and lat_2 are the locations at which
# the cone intersects the earth, so distortion is minimized at those latitudes
alb_70 <- "+proj=aea +lat_1=30 +lat_2=50 +lon_0=-70"
# custom sterographic projection for this nw atlantic, centered using lon_0
stereographic_north <- "+proj=stere +lat_0=90 +lat_ts=75 +lon_0=-57"
# equal earth projection - not working, don't think sf has functionality for it yet
eqearth_proj <- "+proj=eqearth"
# Choose crs using parameter
projection_crs <- switch(
tolower(params$region),
"gulf of maine" = alb_70,
"cpr gulf of maine" = alb_70,
"northwest atlantic" = stereographic_north)
# Transform all the polygons
region_projected <- st_transform(region_extent, crs = projection_crs)
canada_projected <- st_transform(canada, crs = projection_crs)
newengland_projected <- st_transform(new_england, crs = projection_crs)
greenland_projected <- st_transform(greenland, crs = projection_crs)
# coord_sf Crop bounds in projection units for coord_sf
crop_x <- st_bbox(region_projected)[c(1,3)]
crop_y <- st_bbox(region_projected)[c(2,4)]
# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-73185.13, 386673.34)
crop_y <- c(4269044, 4813146)}
# Lower the ymin a touch for NW Atlantic
if(tolower(params$region) == "northwest atlantic"){crop_y <- crop_y - c(100000, 0)}
# # Clip Raster - Convert to stars
# nw_ras <- crop(rotate(july_avg), extent(-73, -40, 40.5, 70))
# nw_st <- st_as_stars(nw_ras)
# Warp to grid of chosen CRS
projection_grid <- region_st %>%
st_transform(projection_crs) %>%
st_bbox() %>%
st_as_stars()
region_warp_ras <- region_st %>%
st_warp(projection_grid)
# Plot everything together
ggplot() +
geom_stars(data = region_warp_ras) +
#geom_sf(data = region_projected, color = "black", fill = "transparent", linetype = 2) +
geom_sf(data = newengland_projected, fill = "gray90") +
geom_sf(data = canada_projected, fill = "gray90") +
geom_sf(data = greenland_projected, fill = "gray90") +
coord_sf(crs = projection_crs,
xlim = crop_x, ylim = crop_y, expand = T) +
map_theme +
scale_fill_gradient2(low = "blue",
mid = "white",
high = "red",
na.value = "transparent") +
guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
theme(
panel.border = element_rect(color = "black", fill = NA),
plot.background = element_rect(color = "transparent", fill = "transparent"),
axis.title.x = element_blank(), # turn off titles
axis.title.y = element_blank(),
legend.position = "bottom",
legend.title.align = 0.5)
Warming Rates
Global warming rates are calculated by getting the average temperature for each year across all cells in the data, then using a linear regression for each to determine the annual warming rate in degrees celsius.
The warming rates are then ranked from high to low to identify areas with the most rapid ocean warming. The following maps display the warming rates and the percentile ranks for the top 20% of warming rates, or rates above the 80th percentile of the data.
# Access information to netcdf on box
rates_path <- shared.path(group = "NSF OKN", folder = "oisst/warming_rates/")
rates_stack <- stack(str_c(rates_path, "annual_warming_rates.nc"), varname = "annual_warming_rate")
ranks_stack <- stack(str_c(rates_path, "annual_warming_rates.nc"), varname = "rate_percentile")
# Going to reclassify raster:
# Source: https://www.earthdatascience.org/courses/earth-analytics/lidar-raster-data-r/classify-raster/
# create classification matrix
reclass_df <- c(0.00, 0.05, 0,
0.05, 0.10, 5,
0.10, 0.15, 10,
0.15, 0.20, 15,
0.20, 0.25, 20,
0.25, 0.30, 25,
0.30, 0.35, 30,
0.35, 0.40, 35,
0.40, 0.45, 40,
0.45, 0.50, 45,
0.50, 0.55, 50,
0.55, 0.60, 55,
0.60, 0.65, 60,
0.65, 0.70, 65,
0.70, 0.75, 70,
0.75, 0.80, 75,
0.80, 0.85, 80,
0.85, 0.90, 85,
0.90, 0.95, 90,
0.95, 1, 95)
# reshape the object into a matrix with columns and rows
reclass_m <- matrix(reclass_df,
ncol = 3,
byrow = TRUE)
# reclassify the raster using the reclass object - reclass_m
ranks_classified <- reclassify(ranks_stack,
reclass_m)
# # Masking Below percentile cutoff - for ranking raster and rates raster
percentile_cutoff <- 80
ranks_stack[ranks_classified < percentile_cutoff] <- NA
rates_stack[ranks_classified < percentile_cutoff] <- NA
# plot(ranks_classified, main = "Reclassified Warming Ranks - top 20%")
# plot(rates_stack, main = "Top 20% Warming Rates")
# Converting to stars
rates_st <- st_as_stars(rotate(rates_stack))
ranks_st <- st_as_stars(rotate(ranks_stack))
ranks_c_st <- st_as_stars(rotate(ranks_classified))
#get scale ranges so they still pop
rate_range <- c("min" = cellStats(rates_stack, 'min')[1], "max" = cellStats(rates_stack,'max')[1])Global Rates
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown."))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown."))
A work by Adam A. Kemberling
Akemberling@gmri.org